fix_results <- do.call('rbind', lapply(list.files(pattern = "demo-sum-fixed"), read.csv))
write.csv(fix_results, "All-fixed.csv", row.names = FALSE)

marg_fixed <- do.call('rbind', lapply(list.files(pattern = "demo-marg-fixed"), read.csv))
write.csv(marg_fixed, "All-marg-fixed.csv", row.names = FALSE)

kappa_results <- do.call('rbind', lapply(list.files(pattern = "demo-marg-kappa"), read.csv))
write.csv(kappa_results, "All-marg-kappa.csv", row.names = FALSE)

range_results <- do.call('rbind', lapply(list.files(pattern = "demo-marg-range"), read.csv))
write.csv(range_results, "All-marg-range.csv", row.names = FALSE)

var_results <- do.call('rbind', lapply(list.files(pattern = "demo-marg-var"), read.csv))
write.csv(var_results, "All-marg-var.csv", row.names = FALSE)

dic_results <- do.call('rbind', lapply(list.files(pattern = "demo-dic"), read.csv))
write.csv(dic_results, "All-DIC.csv", row.names = FALSE)

waic_results <- do.call('rbind', lapply(list.files(pattern = "demo-waic"), read.csv))
write.csv(waic_results, "All-WAIC.csv", row.names = FALSE)

time_results <- do.call('rbind', lapply(list.files(pattern = "demo-time"), read.csv))
write.csv(time_results, "All-times.csv", row.names = FALSE)

Read Data

levels <- paste0("Mesh", 1:36)
fixed <- read.csv("All-fixed.csv")
fixed$meshes <- factor(fixed$Mesh, levels = levels)

margfixed <- read.csv("All-marg-fixed.csv")
margfixed$meshes <- factor(margfixed$Mesh, levels = levels)

margkappa <- read.csv("All-marg-kappa.csv")
margkappa$meshes <- factor(margkappa$Mesh, levels = levels)

margrange <- read.csv("All-marg-range.csv")
margrange$meshes <- factor(margrange$Mesh, levels = levels)

margvar <- read.csv("All-marg-var.csv")
margvar$meshes <- factor(margvar$Mesh, levels = levels)

dic <- read.csv("All-DIC.csv")
dic$meshes <- factor(dic$Mesh, levels = levels)
subset(dic, is.na(dic$DIC))
## [1] X       DIC     Set     Pattern Mesh    meshes 
## <0 rows> (or 0-length row.names)
waic <- read.csv("All-WAIC.csv")
waic$meshes <- factor(waic$Mesh, levels = levels)
subset(waic, is.na(waic$WAIC))
## [1] X       WAIC    Set     Pattern Mesh    meshes 
## <0 rows> (or 0-length row.names)
times <- read.csv("All-times.csv")
times$meshes <- factor(times$Mesh, levels = levels)

Summary fixed

MSE

fixed$diff <- fixed$True - fixed$mean
MSE <- fixed %>%
  group_by(meshes) %>%
  summarize(across("diff", ~ sum(.x^2)/n()))
kable(MSE, align = "cc", digits = 4)
meshes diff
Mesh1 0.1319
Mesh2 0.1338
Mesh3 9.2921
Mesh4 0.1275
Mesh5 0.1196
Mesh6 0.1174
Mesh7 0.1188
Mesh8 0.1306
Mesh9 0.1236
Mesh10 0.1256
Mesh11 0.1544
Mesh12 0.1273
Mesh13 0.1148
Mesh14 0.1130
Mesh15 0.1117
Mesh16 0.1135
Mesh17 0.1048
Mesh18 0.1047
Mesh19 0.1042
Mesh20 0.1040
Mesh21 0.1039
Mesh22 0.1046
Mesh23 0.1042
Mesh24 0.1044
Mesh25 0.0943
Mesh26 0.0936
Mesh27 0.0931
Mesh28 0.0935
Mesh29 0.0942
Mesh30 0.0936
Mesh31 0.0932
Mesh32 0.0937
Mesh33 0.0922
Mesh34 0.0934
Mesh35 0.0933
Mesh36 0.0926

MSE by Set

fixed.set <- split(fixed, fixed$Set)
MSE1 <- fixed.set[[1]] %>%
  group_by(meshes) %>%
  summarize(across("diff", ~ sum(.x^2)/n()))
kable(MSE1, align = "cc", digits = 6)
meshes diff
Mesh1 0.170875
Mesh2 0.174495
Mesh3 11.392104
Mesh4 0.165681
Mesh5 0.140124
Mesh6 0.139521
Mesh7 0.148920
Mesh8 0.177168
Mesh9 0.154710
Mesh10 0.156810
Mesh11 0.213650
Mesh12 0.171066
Mesh13 0.129795
Mesh14 0.127075
Mesh15 0.126636
Mesh16 0.129047
Mesh17 0.103441
Mesh18 0.103829
Mesh19 0.109410
Mesh20 0.106919
Mesh21 0.102797
Mesh22 0.103420
Mesh23 0.109289
Mesh24 0.106897
Mesh25 0.079862
Mesh26 0.078496
Mesh27 0.074512
Mesh28 0.072996
Mesh29 0.079520
Mesh30 0.077955
Mesh31 0.074525
Mesh32 0.072976
Mesh33 0.074474
Mesh34 0.074608
Mesh35 0.073448
Mesh36 0.071083 
MSE1 <- MSE1[-3, ]

MSE2 <- fixed.set[[2]] %>%
  group_by(meshes) %>%
  summarize(across("diff", ~ sum(.x^2)/n()))
kable(MSE2, align = "cc", digits = 6)
meshes diff
Mesh1 0.047330
Mesh2 0.046173
Mesh3 25.437216
Mesh4 0.047109
Mesh5 0.044257
Mesh6 0.042303
Mesh7 0.048292
Mesh8 0.061501
Mesh9 0.042383
Mesh10 0.040082
Mesh11 0.111418
Mesh12 0.051091
Mesh13 0.045199
Mesh14 0.043948
Mesh15 0.051148
Mesh16 0.053905
Mesh17 0.048535
Mesh18 0.046310
Mesh19 0.053251
Mesh20 0.056425
Mesh21 0.048683
Mesh22 0.046425
Mesh23 0.053358
Mesh24 0.056459
Mesh25 0.052803
Mesh26 0.051787
Mesh27 0.062654
Mesh28 0.064313
Mesh29 0.052948
Mesh30 0.051877
Mesh31 0.062688
Mesh32 0.065502
Mesh33 0.057068
Mesh34 0.052484
Mesh35 0.063235
Mesh36 0.065387 
MSE2 <- MSE2[-3, ]

MSE3 <- fixed.set[[3]] %>%
  group_by(meshes) %>%
  summarize(across("diff", ~ sum(.x^2)/n()))
kable(MSE3, align = "cc", digits = 6)
meshes diff
Mesh1 0.242857
Mesh2 0.248523
Mesh3 0.251470
Mesh4 0.241928
Mesh5 0.234481
Mesh6 0.231593
Mesh7 0.230426
Mesh8 0.234787
Mesh9 0.238443
Mesh10 0.245307
Mesh11 0.239114
Mesh12 0.236907
Mesh13 0.229452
Mesh14 0.227072
Mesh15 0.222277
Mesh16 0.224676
Mesh17 0.212736
Mesh18 0.214108
Mesh19 0.206439
Mesh20 0.205514
Mesh21 0.209891
Mesh22 0.213951
Mesh23 0.206251
Mesh24 0.207265
Mesh25 0.188449
Mesh26 0.188050
Mesh27 0.184298
Mesh28 0.183679
Mesh29 0.188487
Mesh30 0.188721
Mesh31 0.184552
Mesh32 0.183575
Mesh33 0.183126
Mesh34 0.185134
Mesh35 0.184615
Mesh36 0.181190 
MSE4 <- fixed.set[[4]] %>%
  group_by(meshes) %>%
  summarize(across("diff", ~ sum(.x^2)/n()))
kable(MSE4, align = "cc", digits = 6)
meshes diff
Mesh1 0.066483
Mesh2 0.065988
Mesh3 0.087592
Mesh4 0.055456
Mesh5 0.059544
Mesh6 0.056285
Mesh7 0.047668
Mesh8 0.048914
Mesh9 0.058765
Mesh10 0.060050
Mesh11 0.053242
Mesh12 0.050172
Mesh13 0.054791
Mesh14 0.053819
Mesh15 0.046909
Mesh16 0.046484
Mesh17 0.054304
Mesh18 0.054702
Mesh19 0.047767
Mesh20 0.047024
Mesh21 0.054222
Mesh22 0.054533
Mesh23 0.047785
Mesh24 0.047122
Mesh25 0.055987
Mesh26 0.056003
Mesh27 0.050940
Mesh28 0.052878
Mesh29 0.055875
Mesh30 0.055990
Mesh31 0.050862
Mesh32 0.052830
Mesh33 0.054060
Mesh34 0.061345
Mesh35 0.051916
Mesh36 0.052688

The mean square error between the estimated and true mean.

for (i in 1:4){
  mse <- ggplot(get(paste0("MSE", i)), aes(x=meshes, y=diff)) +
    geom_point() +
    scale_color_brewer(palette = "Dark2") +
    theme_linedraw() +
    theme(axis.title.x=element_blank(),
          axis.text.x=element_text(angle = 35, hjust = 1)) +
    labs(x="Mesh",y="MSE", title="The mean square error between the estimated and true mean.")
  print(mse)
}

Summary fixed: fixed by meshes

ggplot(fixed, aes(meshes, mean, group = meshes, color = Set)) +
  geom_point(position = position_jitterdodge(jitter.width = .2, dodge.width = .7), 
             alpha = .5) +
  stat_summary(fun = mean, na.rm = TRUE, 
               geom = "point", shape = "diamond",
               size = 2, color = "black", 
               position = position_dodge(width = .7)) +
  stat_summary(fun.data = mean_cl_normal, na.rm = TRUE, 
               geom = "errorbar", width = .2, color = "black",
               position = position_dodge(width = .7)) +
  scale_color_brewer(palette = "Dark2") +
  theme_linedraw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(angle = 35, hjust = 1))

fixed by Set

ggplot(fixed, aes(meshes, mean, group = Set, color = Set)) +
  geom_point(position = position_jitterdodge(jitter.width = .2, dodge.width = .7), 
             alpha = .5) +
  stat_summary(fun = mean, na.rm = TRUE, 
               geom = "point", shape = "diamond",
               size = 2, color = "black", 
               position = position_dodge(width = .7)) +
  stat_summary(fun.data = mean_cl_normal, na.rm = TRUE, 
               geom = "errorbar", width = .2, color = "black",
               position = position_dodge(width = .7)) +
  scale_color_brewer(palette = "Dark2") +
  theme_linedraw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(angle = 35, hjust = 1))

The marginal fixed

The posterior distribution of the log-Cox model parameters.

margfixed.set <- split(margfixed, margfixed$Set)
truefixed <- c(4,4,2,2)
### Set 1
ggplot(margfixed.set[[1]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = truefixed[1], col = 1, linetype = "dashed") +
  labs(x=expression(beta[0]), y="Density", title = paste("Set", 1, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 2
ggplot(margfixed.set[[2]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = truefixed[2], col = 1, linetype = "dashed") +
  labs(x=expression(beta[0]), y="Density", title = paste("Set", 2, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 3
ggplot(margfixed.set[[3]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = truefixed[3], col = 1, linetype = "dashed") +
  labs(x=expression(beta[0]), y="Density", title = paste("Set", 3, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 4
ggplot(margfixed.set[[4]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = truefixed[4], col = 1, linetype = "dashed") +
  labs(x=expression(beta[0]), y="Density", title = paste("Set", 4, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

marginals.kappa

margkappa.set <- split(margkappa, margkappa$Set)
truekappa <- c(10,2,10,2)

### Set 1
ggplot(margkappa.set[[1]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = truekappa[1], col = 1, linetype = "dashed") +
  labs(x=expression(kappa), y="Density", title = paste("Set", 1, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 2
ggplot(margkappa.set[[2]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = truekappa[2], col = 1, linetype = "dashed") +
  labs(x=expression(kappa), y="Density", title = paste("Set", 2, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 3
ggplot(margkappa.set[[3]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = truekappa[3], col = 1, linetype = "dashed") +
  labs(x=expression(kappa), y="Density", title = paste("Set", 3, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 4
ggplot(margkappa.set[[4]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = truekappa[4], col = 1, linetype = "dashed") +
  labs(x=expression(kappa), y="Density", title = paste("Set", 4, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

marginals.variance.nominal

margvar.set <- split(margvar, margvar$Set)
true <- c(1,1,1,1)
### Set 1
ggplot(margvar.set[[1]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = true[1], col = 1, linetype = "dashed") +
  labs(x=expression(sigma^2), y="Density", title = paste("Set", 1, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 2
ggplot(margvar.set[[2]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = true[2], col = 1, linetype = "dashed") +
  labs(x=expression(sigma^2), y="Density", title = paste("Set", 2, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 3
ggplot(margvar.set[[3]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = true[3], col = 1, linetype = "dashed") +
  labs(x=expression(sigma^2), y="Density", title = paste("Set", 3, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 4
ggplot(margvar.set[[4]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = true[4], col = 1, linetype = "dashed") +
  labs(x=expression(sigma^2), y="Density", title = paste("Set", 4, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

marginals.range.nominal

margrange.set <- split(margrange, margrange$Set)
### Set 1
ggplot(margrange.set[[1]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = sqrt(8)/truekappa[1], col = 1, linetype = "dashed") +
  labs(x="Nominal range", y="Density", title = paste("Set", 1, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 2
ggplot(margrange.set[[2]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = sqrt(8)/truekappa[2], col = 1, linetype = "dashed") +
  labs(x="Nominal range", y="Density", title = paste("Set", 2, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 3
ggplot(margrange.set[[3]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = sqrt(8)/truekappa[3], col = 1, linetype = "dashed") +
  labs(x="Nominal range", y="Density", title = paste("Set", 3, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

### Set 4
ggplot(margrange.set[[4]], aes(x=x, y=y, group=Pattern, color=Pattern)) + 
  geom_line() + 
  facet_wrap(~ meshes, scales = "free") +
  geom_vline(xintercept = sqrt(8)/truekappa[4], col = 1, linetype = "dashed") +
  labs(x="Nominal range", y="Density", title = paste("Set", 4, sep=" ")) +
  scale_color_brewer(palette = "Set1") +
  theme_light() +
  theme(plot.title = element_text(size=11),
        axis.text.x = element_text(hjust = 1, size = unit(7, "lines")),
        axis.text.y = element_text(hjust = 1, size = unit(7, "lines")),
        panel.spacing = unit(0.5, "lines"), 
        strip.text = element_text(size=8, face="bold"))

DIC

DIC by meshes

dic_set <- split(dic, dic$Set)
dic_set1 <- dic_set[[1]] %>%
  group_by(meshes) %>%
  summarise_at(vars(DIC), list(name = mean))
dic_set1[which.min(dic_set1$name),]
## # A tibble: 1 x 2
##   meshes   name
##   <fct>   <dbl>
## 1 Mesh3  -6031.
dic_set1[which.max(dic_set1$name),]
## # A tibble: 1 x 2
##   meshes   name
##   <fct>   <dbl>
## 1 Mesh12 -5589.
dic_set2 <- dic_set[[2]] %>%
  group_by(meshes) %>%
  summarise_at(vars(DIC), list(name = mean))
dic_set2[which.min(dic_set2$name),]
## # A tibble: 1 x 2
##   meshes   name
##   <fct>   <dbl>
## 1 Mesh3  -6880.
dic_set2[which.max(dic_set2$name),]
## # A tibble: 1 x 2
##   meshes   name
##   <fct>   <dbl>
## 1 Mesh32 -6177.
dic_set3 <- dic_set[[3]] %>%
  group_by(meshes) %>%
  summarise_at(vars(DIC), list(name = mean))
dic_set3[which.min(dic_set3$name),]
## # A tibble: 1 x 2
##   meshes  name
##   <fct>  <dbl>
## 1 Mesh34 -335.
dic_set3[which.max(dic_set3$name),]
## # A tibble: 1 x 2
##   meshes  name
##   <fct>  <dbl>
## 1 Mesh2  -330.
dic_set4 <- dic_set[[4]] %>%
  group_by(meshes) %>%
  summarise_at(vars(DIC), list(name = mean))
dic_set4[which.min(dic_set4$name),]
## # A tibble: 1 x 2
##   meshes  name
##   <fct>  <dbl>
## 1 Mesh32 -403.
dic_set4[which.max(dic_set4$name),]
## # A tibble: 1 x 2
##   meshes  name
##   <fct>  <dbl>
## 1 Mesh1  -375.
ggplot(dic, aes(x=meshes, y=DIC, group = meshes, color = Set)) +
  geom_point(position = position_jitterdodge(jitter.width = .2, dodge.width = .7), 
             alpha = .5) +
  stat_summary(fun = mean, na.rm = TRUE, 
               geom = "point", shape = "diamond",
               size = 2, color = "black", 
               position = position_dodge(width = .7)) +
  stat_summary(fun.data = mean_cl_normal, na.rm = TRUE, 
               geom = "errorbar", width = .2, color = "black",
               position = position_dodge(width = .7)) +
  scale_color_brewer(palette = "Dark2") +
  theme_linedraw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(angle = 45, hjust = 1))

DIC by Set

ggplot(dic, aes(x=meshes, y=-DIC, group = Set, color = Set)) +
  geom_point(position = position_jitterdodge(jitter.width = .2, dodge.width = .7), 
             alpha = .5) +
  stat_summary(fun = mean, na.rm = TRUE, 
               geom = "point", shape = "diamond",
               size = 2, color = "black", 
               position = position_dodge(width = .7)) +
  stat_summary(fun.data = mean_cl_normal, na.rm = TRUE, 
               geom = "errorbar", width = .2, color = "black",
               position = position_dodge(width = .7)) +
  scale_color_brewer(palette = "Dark2") +
  theme_linedraw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(angle = 45, hjust = 1))

WAIC

WAIC by meshes

waic_set <- split(waic, waic$Set)
waic_set[[1]][which.min(-waic_set[[1]]$WAIC),]
##    X      WAIC  Set  Pattern   Mesh meshes
## 42 2 -4419.122 Set1 Pattern2 Mesh11 Mesh11
waic_set[[2]][which.min(-waic_set[[2]]$WAIC),]
##     X      WAIC  Set  Pattern  Mesh meshes
## 446 6 -2454.885 Set2 Pattern6 Mesh3  Mesh3
waic_set[[3]][which.min(-waic_set[[3]]$WAIC),]
##     X      WAIC  Set   Pattern  Mesh meshes
## 11 11 -250.8282 Set3 Pattern11 Mesh1  Mesh1
waic_set[[4]][which.min(-waic_set[[4]]$WAIC),]
##      X    WAIC  Set   Pattern  Mesh meshes
## 456 16 -98.049 Set4 Pattern16 Mesh3  Mesh3
ggplot(waic, aes(x=meshes, y=-WAIC, group = meshes, color = Set)) +
  geom_point(position = position_jitterdodge(jitter.width = .2, dodge.width = .7), 
             alpha = .5) +
  stat_summary(fun = mean, na.rm = TRUE, 
               geom = "point", shape = "diamond",
               size = 2, color = "black", 
               position = position_dodge(width = .7)) +
  stat_summary(fun.data = mean_cl_normal, na.rm = TRUE, 
               geom = "errorbar", width = .2, color = "black",
               position = position_dodge(width = .7)) +
  scale_color_brewer(palette = "Dark2") +
  theme_linedraw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(angle = 45, hjust = 1))

WAIC by Set

ggplot(waic, aes(x=meshes, y=-WAIC, group = Set, color = Set)) +
  geom_point(position = position_jitterdodge(jitter.width = .2, dodge.width = .7), 
             alpha = .5) +
  stat_summary(fun = mean, na.rm = TRUE, 
               geom = "point", shape = "diamond",
               size = 2, color = "black", 
               position = position_dodge(width = .7)) +
  stat_summary(fun.data = mean_cl_normal, na.rm = TRUE, 
               geom = "errorbar", width = .2, color = "black",
               position = position_dodge(width = .7)) +
  scale_color_brewer(palette = "Dark2") +
  theme_linedraw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(angle = 45, hjust = 1))

Computational time

By Mesh

# stat_summary with arg "fun.data": 
# A function that is given the complete data and should return a data frame
# with variables ymin, y, and ymax (for use in plotting ranges).
# mean_cl_normal() is intended for use with stat_summary. It calculates
# sample mean and lower and upper Gaussian confidence limits based on the t-distribution
ggplot(times, aes(x=meshes, y=elapsed, group=meshes, color=Set)) +
  geom_point(position = position_jitterdodge(jitter.width = .2, dodge.width = .7), 
             alpha = .5) +
  stat_summary(fun = mean, na.rm = TRUE, 
               geom = "point", color = "black", 
               size = 2, shape = "diamond") +
  stat_summary(fun.data = mean_cl_normal, na.rm =TRUE, 
               geom = "errorbar", width = .2, color = "black") +
  theme(axis.text.x=element_text(angle = 90, size = 7)) +
  scale_color_brewer(palette = "Dark2") +
  theme_linedraw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(angle = 35, hjust = 1))

By Set

ggplot(times, aes(x=meshes, y=elapsed, group=Set, color=Set)) +
  geom_point(position = position_jitterdodge(jitter.width = .2, dodge.width = .7), 
             alpha = .5) +
  stat_summary(fun = mean, na.rm = TRUE, 
               position = position_dodge(width = .7),
               geom = "point", color = "black", 
               size = 2, shape = "diamond") +
  stat_summary(fun.data = mean_cl_normal, na.rm =TRUE, 
               position = position_dodge(width = .7),
               geom = "errorbar", width = .2, color = "black") +
  theme(axis.text.x=element_text(angle = 90, size = 7)) +
  scale_color_brewer(palette = "Dark2") +
  theme_linedraw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(angle = 35, hjust = 1))